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1.  Introduction 
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In  this  paper  we  investigate  the  problem  of  passive  navigation  using  optical  flow 
information.  Suppose  we  are  viewing  a  film.  We  wish  to  determine  the  motion  of  the 
camera  from  the  sequence  of  images,  assuming  that  the  instantaneous  velocity  of  the 
brightness  patterns,  also  called  the  optical  flow,  is  known  at  each  point  in  the  image. 
Several  schemata  for  computing  optical  flow  have  been  suggested  (e.g.  [2],  [3],  [5]). 
Other  papers  (e.g.  [9],  [10],  [11])  have  previously  addressed  the  problem  of  passive 
navigation.  Three  approaches  can  be  taken  towards  a  solution  which  we  term  the 
discrete,  the  differential  and  the  continuous  approach. 

In  the  discrete  approach,  information  about  the  movement  of  brightness  patterns 
at  only  a  few  points  is  used  to  determine  the  motion  of  the  camera.  In  particular,  using 
such  an  approach,  one  attempts  to  identify  and  match  discrete  points  in  a  sequence  of 
images.  Of  interest  in  this  case  is  the  photogrammetric  problem  of  determining  what 
the  minimum  number  of  points  is  from  which  the  motion  can  be  calculated  for  a  given 
number  of  images  [10],  [11],  [14].  This  approach  requires  that  one  tracks  features,  or 
identifies  corresponding  features  in  images  taken  at  different  times. 

In  the  differential  approach,  the  first  and  second  spatial  partial  derivatives  of  the 
optical  flow  are  used  to  compute  the  motion  of  a  camera  [6],  [9].  It  has  been  claimed 
that  it  is  sufficient  to  know  the  optical  flow  and  both  its  first  and  second  derivatives 
at  a  single  point  to  uniquely  determine  the  motion  [9].  This  is  incorrect  (except  for 
a  special  case)  [1].  Furthermore,  noise  in  the  measured  optical  flow  is  accentuated  by 
differentiation. 

In  the  continuous  approach,  the  whole  optical  flow  field  is  used.  A  major  shortcom¬ 
ing  of  both  the  local  and  differential  approaches  is  that  neither  allows  for  errors  in  the 
optical  flow  data.  This  is  why  we  choose  the  continuous  approach  and  devise  a  least- 
squares  technique  to  determine  the  motion  of  the  camera  from  the  measured  optical 
flow.  The  proposed  algorithm  takes  the  abundance  of  available  data  into  account  and 
is  robust  enough  to  allow  numerical  implementation. 

2.  Technical  Prerequisites 

In  this  section  we  review  the  equations  describing  the  relation  between  the  motion 
of  a  camera  and  the  optical  flow  generated.  We  use  essentially  the  same  notation  as 
[9].  A  camera  is  assumed  to  move  through  a  static  environment.  Let  a  coordinate 
system  X ,  Y ,  Z  be  fixed  with  respect  to  the  camera,  with  the  Z-axis  pointing  along 
the  optical  axis.  If  we  wish,  we  can  think  of  the  environment  as  moving  in  relation 
to  this  coordinate  system.  Any  rigid  body  motion  can  be  resolved  into  two  factors, 
a  translation  and  a  rotation.  We  will  denote  by  f  the  translational  component  of 
the  motion  of  the  camera  and  by  Q  its  angular  velocity  (see  also  Figure  1  which  is 
redrawn  from  [9]).  Let  the  instantaneous  coordinates  of  a  point  P  in  the  environment 

be  (X,r,Z). 
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(Note  that  Z  >  0  for  points  in  front  of  the  imaging  system.)  Let  f  be  the  vector 
(XtYfZ)Tf  where  r  denotes  the  transpose  of  a  vector,  then  the  velocity  of  P  with 
respect  to  the  X,Y,  Z  coordinate  system,  is: 

V  =  -T-axf.  (l) 

We  define  the  components  of  t  and  Q  as: 

T  =  (U,V,W)T  G>  =  (A,  B,  C)T .  (2) 

Thus  we  can  rewrite  (1)  in  component  form: 

X'  =  -U  -  BZ  +  CY 

Y'  =  —V  —  CX  +  AZ  (3) 

Z*  =  -W  -AY  +  BX. 

where  9  denotes  differentiation  with  respect  to  time. 

The  optical  Bow  at  each  point  in  the  image  plane  is  the  instantaneous  velocity  of 
the  brightness  pattern  at  that  point.  Let  (x,y)  denote  the  coordinates  of  a  point  in  the 
image  plane  (see  Figure  1).  Since  we  assume  perspective  projection  between  an  object 
point  P  and  the  corresponding  image  point  p,  the  coordinates  of  p  are: 
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y=z- 


The  optical  flow,  denoted  by  (u,  v),  at  a  point  ( x ,  y)  is: 

u  =  x'  v  =  yl.  (5) 

Differentiating  (4)  with  respect  to  time  and  using  (3)  we  obtain  the  following  equations 
for  the  optical  flow: 

V'  Y  7*  TT  W 

U=ZY~  =  -B  +  Cy)- x{--£  -Ay  +  Bx)  (6) 

Y<  YZ'  V  W 

v  =  —  -—T  =[---Cx  +  A)-y{— --Ay  +  Bx). 


We  can  write  these  equations  in  the  form: 

u  =  u,  +  ur 


V  =  Vt  -f  Vr 


where  (ut,vt)  denotes  the  translational  component  of  the  optical  flow  and  (ur,vr)  the 
rotational  component: 

-U  +  xW  -V  +  yW 

«t  = - - -  vt  =  - -  (8) 

ur  =  Axy  —  B(x 2  +  1)  +  Cy  vr  ~  A(y 2  +  1)  —  Bxy  —  Cx . 

So  far  we  have  considered  a  single  point  P.  To  define  the  optical  flow  globally  we 
assume  that  P  lies  on  a  surface  defined  by  a  function  Z  =  Z(X,Y)  which  is  positive  for 
all  values  of  X  and  Y .  With  any  surface  and  any  motion  of  a  camera  we  can  therefore 
associate  a  certain  optical  flow  and  we  say  that  the  surface  and  the  motion  generate 
this  optical  flow. 

Optical  flow,  therefore,  depends  upon  the  six  parameters  of  motion  of  the  camera 
and  upon  the  surface  whose  images  are  analyzed.  Can  all  these  unknowns  be  uniquely 
recaptured  solely  from  optical  flow?  The  answer  is  no.  To  see  this,  consider  a  surface 
S2  which  is  a  dilation  by  a  factor  k  of  a  surface  S\.  Further,  let  two  motions  denoted  by 
M\  and  M2  have  the  same  rotational  component  and  let  their  translational  components 
be  proportional  to  each  other  by  the  same  factor  k  (we  will  say  that  M\  and  M2  are 
similar).  Then  the  optical  flow  generated  by  Si  and  Mi  is  the  same  as  the  optical  flow 
generated  by  S2  and  M%.  This  follows  directly  from  the  definition  of  optical  flow  (8). 
It  is  still  an  open  question  whether  there  are  any  other  pairs  of  distinct  surfaces  and 
motions  which  generate  the  same  optical  flow. 

Determining  the  motion  of  a  camera  from  optical  flow  is  much  easier  if  we  are  told 
that  the  motion  is  purely  translational  or  purely  rotational.  In  the  next  two  sections 
we  will  deal  with  these  two  special  cases.  Then  w c  shall  analyze  the  case  where  no  a 
priori  assumptions  about  the  motion  of  the  camera  are  made. 


3.  Translational  Case 


5 


In  this  section  we  discuss  the  case  where  the  motion  of  the  camera  is  assumed  to 
be  purely  translational.  As  before,  let  f  =  (U ,  V,  IV)  be  the  velocity  of  the  camera. 
Then  the  following  equations  hold  (see  (8)): 


-U  +  xW 


3.1.  Similar  Surfaces  and  Similar  Motions 


-V+y  W 


It  will  be  shown  next  that  if  two  flows  generate  the  same  optical  flow,  and  we  know 
that  the  motions  are  purely  translational,  then  the  two  surfaces  are  similar  and  the  two 
camera  motions  are  similar.  Let  Z\  and  Zi  be  two  surfaces  and  let  7\  =  [U\ ,  Vi,  Wi)T 
and  7j  —  (U2,  Vj,  W2)T  define  two  different  motions  of  a  camera,  such  that  Z\  and  T2 
and  Z2  and  7s  generate  the  same  optical  flow: 


-Ui+*WX 

Zi 

-U2  +  xW3 

z2 


-vi-fyiv! 

Z\  1 
-Va  +  yWa 
Z2 


Eliminating  Zlt  Z2,  u  and  v  from  these  equations  we  obtain: 

— U\  +  xW2  -U2  +  xW2 
-Vl+yWl  ~  - V2  -f  yW2 

We  can  rewrite  this  equation  as: 

(~Ux  +  xW1)(-V2  +  yW2)  =  (~U2  +  xWjK-Vj  +  yWt), 


UiV- 2  -  xV2W\  -  yUiW2  +  xyWiWa  =  U2V ,  -  xViW2  -  yU2W2  +  xyW2Wx.  (14) 

Since  we  assumed  that  Zj  and  7j  and  Z2  and  t2  generate  the  same  optical  flow,  the 
above  equation  must  hold  for  ail  x  and  y.  Therefore  the  following  equations  have  to 
hold: 

Ulv2  =  u2vt 

-v2  Wy  =  -VxW2  (15) 

-tWs  =  -U2WX. 

These  equations  can  be  rewritten  as: 

Uf.Vx-.Wt  =  U2:V2:W2  (16) 


m&zr  > 
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from  which  it  follows  that  Z 2  is  a  dilation  of  Zx.  It  is  clear  that  the  scaling  factor 
between  Z\  and  Z2  (or  equivalently  between  Tx  and  f 2)  cannot  be  recovered  from  the 
optical  flow,  regardless  of  the  number  of  points  at  which  the  flow  is  known.  By  uniqueiy 
determining  the  motion  of  the  camera,  we  will  mean  that  the  motion  is  uniquely 
determined  up  to  a  constant  scaling  factor. 


3.2.  Least-Squares  Formulation 

In  general,  the  direction  of  the  optical  flow  at  two  points  in  the  image  plane 
determine  the  motion  of  a  camera  in  pure  translation  uniquely.  There  is  a  drawback 
however  to  utilizing  so  little  of  the  available  information.  The  optical  flow  we  measure 
is  corrupted  by  noise  and  it  is  desirable  to  develop  a  robust  method  which  takes  this 
into  account.  Thus  we  suggest  using  a  least-* squares  method  [4],  [12]  to  determine  the 
movement  parameters  and  the  surface  (i.e.,  the  best  fit  with  respect  to  some  norm). 

For  the  following  we  assume  that  the  image  plane  is  the  rectangle  ze[ — w,  in]  and 
y[ — h,h].  The  same  method  applies  if  the  image  has  some  other  shape.  (In  fact,  it 
can  be  used  on  sub-images  corresponding  to  individual  objects  in  the  case  that  the 
environment  contains  objects  which  may  move  relative  to  one  another).  Furthermore 
we  have  to  assume  that  l/Z  is  a  bounded  function  and  that  the  set  of  points  where  1  /Z 
is  discontinuous  is  of  measure  zero.  This  condition  on  l/Z  assures  us  that  all  necessary 
integrations  can  be  carried  out.  We  wish  to  minimize  the  following  expression: 


-U  +  xW  2 

z  1 


+  (t>  — 


- V  +  yW 
Z 


)2)dxdy. 


(17) 


In  this  case  then,  we  determine  the  best  fit  with  respect  to  the  ML 2  norm  which  is 
defined  as: 


II  /(*»  y)  11=  [  (  \f{x,y)]3dxdy. 

J  — hJ  —  W 


(18) 


The  steps  in  the  least-squares  method  are  as  follows:  First  we  determine  that  Z  which 
minimizes  the  integrand  of  (17)  at  every  point  (x,  y).  Then  we  determine  the  values  of 
U,  V  and  W  which  minimize  the  integral  (17). 

Let  us  introduce  the  following  abbreviations: 


a  =  —U  -f  xW  p  =  +  yW. 


Note  that  the  expected  flow,  given  U,  V  and  W  is  simply: 


Then  we  can  rewrite  (17)  as: 


f)2  +  (*■ -  f  )*]**• 


(19) 


(20) 


(21) 


WCT<rilr»)V 


<ua,va) 


L:  u/3  -  v  i  =  0 


Figure  2.  Geometrical  Interpretation 


We  proceed  now  with  the  first  step  of  our  minimization  method.  Differentiating  the 
integrand  of  (17)  with  respect  to  Z  and  setting  the  resulting  expression  equal  to  zero 
yields: 

(u _  +  (v  —  -)A;  =  o.  (22) 

'  Z}Z2^K  Z'Za  V 


v-  zjZ2  '  v  z'Z2 

Therefore  we  can  write  Z  as:  „ 

Z  =  2L±£_.  (23) 

va  -f  VP 

This  equation,  by  the  way,  imposes  a  constraint  on  U,  V  and  W,  since  Z  must  be 
positive.  We  do  not  make  use  of  this  except  to  help  us  pick  amongst  two  opposite 
solutions  for  the  translational  velocity  later  on.  Note  that  now: 


a  -u/3  —  va 
U~  Z~  ^  ot2  +p 

and  we  can  therefore  rewrite  (17)  as: 


j8  u0  —  vo 

v---  -«q2 


/:/: 


(u/9  —  vo)a 


J-J-u,  o’+/J»  ~  ' 

It  should  be  clear,  by  the  way,  that  uniformly  scaling  U,  V  and  W  does  not  change 
the  value  of  (25).  This  is  a  reflection  of  the  fact  that  we  can  determine  the  motion 
parameters  only  up  to  a  constant  factor. 

Before  proceeding  with  the  second  step,  we  give  a  geometrical  interpretation  in 
Figure  2  of  what  we  have  so  far.  Suppose  that  the  motion  parameters  U}  V,  and 


fesSfte**  ■  *  - 
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W  are  given.  At  any  given  point,  say  (xq,  j/o),  optical  flow  depends  not  only  upon  the 
motion  parameters  but  also  upon  the  value  of  Z  at  that  point,  2T0  say.  However,  the 
direction  of  (u,v)  does  not  depend  upon  Zo-  The  point  (u,v)  must  lie  along  the  line  L 
in  the  uv-plane  defined  by  the  equation  ti0  —  va  —  0.  Let  the  measured  optical  flow 
at  (zo.yo)  be  denoted  by  (um,vm),  and  let  the  closest  point  on  the  line  L  be  (ua,va). 
This  corresponds  to  a  particular  Za  given  by  (23).  The  remaining  error  is  the  distance 
between  the  point  (um,  vm)  and  the  line  L.  The  square  of  this  distance  is  given  by  the 
integrand  of  (25). 

For  the  second  step,  we  differentiate  (25)  with  Tespect  to  U,V  and  W  and  Bet  the 
resulting  expressions  equal  to  zero: 


dxdy  —  0 


C£Jv-V+^',a+^=o- 


Let  us  introduce  the  following  abbreviation: 


_  (u/?  —  «a)(tttt  V0) 
(<*2  +  /?2)2  ' 


Then  equations  (26)  can  be  rewritten  as: 


[  [  [(— V  +  yW)K\dxdy  =  0 

—  f  I  ((-If  +  xW)K  ]dxdy  =  0 
J — hJ — w 

I  I  \(~VU  +  xV)K]dxdy  =  0. 
J  —KJ  —W 


The  sum  of  U  times  the  first  integral,  V  times  the  second  integral,  and  W  times  the 
third  integral  is  identically  zero.  Thus  the  three  equations  are  linearly  dependent.  This 
is  to  be  expected,  for  if: 

f(kU,  kV,  kW)  =  f{V,V,  W),  (29) 

where  /  is  a  differentiable  function  and  k  a  constant,  then: 

uW  +  v^  +  wm-°-  <*» 


(30) 


o 
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The  result  is  also  consistent  with  the  fact  that  only  two  equations  are  needed,  since  the 
translational  velocity  can  be  determined  only  up  to  a  constant  factor.  Unfortunately 
equations  (28)  are  nonlinear  in  Ut  V  and  W  and  we  are  not  able  to  show  that  they  have 
a  unique  (up  to  a  constant  scaling  factor)  solution. 

3.3.  Using  a  Different  Norm 

There  is  a  way,  however,  to  devise  a  least-squares  method  which  allows  us  to  display 
a  closed  form  solution  for  the  motion  parameters.  Instead  of  minimizing  (17),  we  will 
try  to  minimize  the  following  expression: 

LL  + <•  -  =^)!K»! + ***  on 

obtained  by  multiplying  the  integrand  of  (17)  by  a2  -f-  /?2-  Then  we  apply  the  same 
least- squares  method  as  before  to  (31).  When  the  measured  optical  flow  is  not  corrupted 
by  noise,  both  (31)  and  (17)  can  be  made  equal  to  zero  by  substituting  the  correct  motion 
parameters.  We  thus  obtain  the  same  solution  for  the  motion  parameters  whether  we 
minimize  (31)  or  (17).  If  the  measured  optical  flow  is  not  exact,  then  using  expression 
(31)  for  our  minimization,  we  obtain  the  best  fit  with  respect  not  to  the  MLa  norm, 
but  to  another  norm  which  we  call  the  MLap  norm: 

II  /(*,»)  Ha*=  /  f  [/(*,  y)]2(“2  +  P3)dxdy.  (32) 

What  we  have  here  is  a  minimization  in  which  the  error  contributions  are  weighted , 
greater  importance  being  given  to  points  where  the  optical  flow  velocity  is  larger.  This 
is  most  appropriate  when  the  measurement  of  larger  velocities  is  more  accurate. 

Which  norm  gives  the  best  results  depends  on  the  properties  of  the  noise  in  the 
measured  optical  flow.  The  first  norm  is  better  suited  to  the  sitation  where  the  noise  in 
the  measurements  is  independent  of  the  magnitude  of  the  optical  flow.  Note  also  that 
if  we  really  want  the  minimum  with  respect  to  the  ML 2  norm,  we  can  use  the  results 
of  the  minimization  with  respect  to  the  MLa$  norm  as  starting  values  in  a  numerical 
minimization. 

We  discuss  now  our  least-squaie^  uiethod  in  the  case  where  the  norm  is  chosen  to 
be  MLa0.  First  we  determine  z,  by  differentiating  the  integrand  of  (31)  with  respect 
to  Z  and  setting  the  result  equal  to  zero.  We  again  get  (22): 

—  z^Z *  —  ~zi~zi  =  °’  ^ 


from  which  it  follows  that  (23): 


a2  ±JP_ 

ua  +  v0' 


(34) 


So  we  want  to  minimize: 
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f  j  (u/3  —  va)2dxdy. 

J  —  hJ  —  ty 

Let  us  call  this  integral  g(U,  V,  W),  then,  since: 

u0  —  va  =  ( vU  —  uV-)  —  (xv  —  yu)W, 


(35) 


(36) 


we  have: 

g{U,  V,  W)  =  aU 2  +  bV2  +  cW 2  +  2 dUV  +  2 eVW  +  2fWU,  (37) 


where: 


fh  rw 

a  —  1  1  v2dxdy 

J  —  kJ—w 

fh  rw 

b  =  1  1  u2dxdy 

rh  /*W 

c  =  /  /  ( xv  —  yxi)2dxdy 

fh  rw 

d  =  —  /  /  uvdxdy 

J  —  hJ  —w 

(38) 

e  =  /  /  u(xv  —  yu)dxdy 

J  —  hJ  —  Ul 

fh  rw 

f  =  —  I  I  v(xv  —  yu)dxdy. 

J — hJ — ty 

Now  g(U ,  V}  W)  cannot  be  negative,  and  g(U,  V,  W)  =  0  for  U  =  V  =  W  =  0.  Thus 
a  minimum  can  be  found  by  inspection,  but  is  not  what  we  might  have  hoped  for.  In 
fact,  to  determine  the  translational  velocity  using  our  least-squares  method  we  have  to 
solve  the  following  homogeneous  equation  for  T : 

O 

II 

(39) 

where  G  is  the  matrix: 

/a  d  f\ 

G  = 

■G : :) 

(40) 

Clearly  (39)  has  a  solution  other  then  2ero  if  and  only  if  the  determinant  of  G  is  zero. 
Then  the  three  equations  (39)  are  linearly  dependent  and  T  can  be  determined  up  to 
a  constant  factor.  In  general,  however,  as  the  data  is  corrupted  by  noise,  g  cannot  be 
made  equal  to  zero  for  non-zero  translational  velocity  and  so  T  =  (0,0, 0)T  will  be  the 
only  solution  to  (39).  To  see  this  in  another  way,  note  that  g  has  the  following  form: 

g{kU9  kV9  kW )  =  k2g(U,  Vt  W)  (41) 

where  k  is  a  constant.  Clearly  g{U,V}W)  assumes  its  minimal  value  for  {7  =  V  = 
W  =  0. 

What  we  are  really  interested  in,  is  determining  the  direction  of  f  which  minimizes 
gf  for  a  fixed  length  of  T.  Hence  we  impose  the  constraint  that  f  be  a  unit  vector. 
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If  T  is  constrained  to  have  unit  magnitude,  the  minimum  value  of  g  is  the  smallest 
eigenvalue  of  the  matrix  G  and  the  value  of  T  for  which  g  assumes  its  minimum  can  be 
found  by  determining  the  eigenvector  corresponding  to  this  eigenvalue  [8].  This  follows 
from  the  observation  that  g  is  a  quadratic  form  which  can  be  written  as: 

g{U,V,W)  =  fTGt.  (42) 

Note  that  G  is  a  positive  semidefinite  hermitian  matrix  as  a  >  0,  6  >  0,  c  >  0,  a6  > 
d2,  6c  >  e2  and  ca  >  /2.  (The  last  three  inequalities  follow  from  the  Cauchy- Schwarz 
inequality  [7],  [8]).  Hence  all  eigenvalues  are  real  and  non-negative  and  are  the  solutions 
X  of  the  third  degree  polynomial: 

X3  —  (a  +  6  +  +  (<*6  +  6c  +  ca  —  d2  —  e2  —  /2)X  +  (43) 

(ae2  +  6/2  +  cd2  -  a6c  —  2<Uf)  =  0. 

There  is  an  explicit  formula  for  the  least  positive  root  in  terms  of  the  real  and  imaginary 
parts  of  the  roots  of  the  quadratic  resolvent  of  the  cubic.  In  our  case  this  gives  us  the 
desired  smallest  root,  since  the  roots  cannot  be  negative.  For  the  sake  of  completeness, 
however,  various  pathological  cases  that  might  come  up  will  be  discussed  next,  even 
though  they  are  of  little  practical  interest. 

Note  that  X  =  0  is  an  eigenvalue  if  and  only  if  G  is  singular,  that  is,  if  the  constant 
term  in  the  polynomial  (43)  equals  zero.  In  fact,  if  the  determinant  of  G  is  zero  one 
can  find  a  velocity  T  which  makes  g  zero.  It  follows  from  a  theorem  in  calculus  that 
this  happens  only  when  the  optical  Bow  is  either  not  corrupted  by  noise  at  all  or  only 
at  a  few  points.  The  theorem  states  that  if  the  integral  of  the  square  of  a  bounded  and 
continuous  function  is  zero  then  the  function  itself  is  zero.  Hence  errors  can  only  occur 
at  points  where  the  optical  flow  is  discontinuous,  and  these  are  exactly  the  points  where 
the  surface  defined  by  Z  is  discontinuous.  (These  are  also  the  places  where  existing 
methods  for  computing  the  optical  flow  [5]  are  subject  to  large  errors). 

It  is  impossible  for  exactly  two  eigenvalues  to  be  zero,  since  this  would  imply  that 
the  coefficient  of  X  in  the  polynomial  (43)  equalled  zero,  while  the  coefficient  of  X2  did 
not.  That  in  turn  would  imply  that  ab  =  d2y  be  =  e2,  and  ca  =  /2,  while  a,  6,  and 
c  are  not  all  zero.  For  equality  to  hold  in  the  Cauchy- Schwarz  inequalities,  however, 
tt  and  v  must  both  be  proportional  to  xv  —  yu.  This  can  only  be  true  (for  all  z  and 
y  in  the  image)  if  u=u=0.  But  then  all  six  integrals  become  zero  and  consequently 
all  three  eigenvalues  are  zero.  This  situation  i3  of  little  interest,  since  it  occurs  only 
when  the  optical  flow  data  is  zero  everywhere.  Then  the  velocity  is  zero  too.  Once 
the  smallest  eigenvalue  is  known,  it  is  straightforward  to  find  the  translational  velocity 
which  best  matches  the  given  data.  To  determine  the  eigenvector  corresponding  to  an 
eigenvalue,  say  Xi,  we  have  to  solve  the  following  set  of  linear  equations: 


12 


{a-\x)U  +  dV  +  fW  =  Q 

dU  +  (b-\!)V  +  eW  =  Q  (44) 

fU  +  tV  +  (c  -  \X)W  =  0. 

As  Xi  is  an  eigenvalue,  equations  (44)  are  linearly  dependent.  Let  us  for  a  moment 
assume  that  all  eigenvalues  are  distinct,  that  is,  the  rank  of  the  matrix  (G —  X/),  where 
I  is  the  identity  matrix,  is  two.  Then  we  can  use  any  pair  of  them  to  solve  for  U,V  in 
terms  of  W  say.  There  are  three  ways  of  doing  this.  Fot  numerical  accuracy  we  may 
add  the  three  results  to  get  the  symmetrical  forms: 

U  =  {b  —  Xx)(c  —  Xi)  —  f{b  —  \x)  —  d(c  —  Xx)  +  e(f  +  d  —  e) 

V  =  (c  —  Xi)(a  -  Xj)  -  d(c  -  Xi)  -  e(a  -  X,)  +  f(d  +  «  -  /)  (45) 

W  =  (a  —  Xi)(6  —  Xi)  —  e(a  —  Xi)  —  f(b  —  Xi)  +  d(e  +  /  —  d). 

Note  that  Xi  will  be  very  small,  if  the  data  is  good,  and  one  may  wish  to  just 

approximate  the  exact  solution  by  using  the  above  equations  with  Xi  set  to  zero.  (Then 
there  is  no  need  to  find  the  eigenvalue).  In  any  case,  the  resulting  velocity  may  now  be 
normalized  so  that  its  magnitude  equals  one.  There  is  one  remaining  difficulty,  arising 
from  the  fact  that  if  T  is  a  solution  to  our  minimization  problem,  so  is  —  T.  Only  one 
of  these  solution  will  correspond  to  positive  values  of  Z  in  equation  (34)  however.  ThiB 
can  be  easily  seen  by  evaluating  (34)  at  some  point  in  the  image.  The  case  where  the 
two  smallest  eigenvalues  are  the  same  will  be  discussed  in  one  of  the  next  paragraphs. 

There  is  a  simple  geometrical  interpretation  of  what  we  have  done  so  far.  To  this 
end  we  consider  the  surface  defined  by  g(U,  V ,  VV|  =  k  where  k  is  a  constant.  Note 
that  we  can  always  find  a  new  coordinate  system  UfV,W  in  which  g(U,  V,  W)  can  be 
written  as: 

Xit/2  +  X2V  +  X3W  =*  (46) 

where  \  for  i  =  1, 2, 3  are  the  three  eigenvalues  of  the  quadratic  form.  If  the  eigenvalues 
are  all  non-zero,  the  surface  g(Ut  V,  W)  =  k  is  an  ellipsoid  with  three  orthogonal  semi¬ 
axes  of  length  y/k/\i.  We  are  particularly  interested  in  the  case  where  the  constant 
k  is  the  smallest  eigenvalue.  Then  all  three  semi-axes  have  lengths  less  than  or  equal 
to  one.  Hence  the  ellipsoid  lies  within  the  unit  sphere.  If  the  two  smallest  eigenvalues 
are  distinct,  the  unit  sphere  touches  the  ellipsoid  in  two  places,  corresponding  to  the 
largest  axis.  If  the  two  smaller  eigenvalues  happen  to  be  the  same,  however,  the  unit 
sphere  touches  the  ellipsoid  along  a  circle  and  as  a  result  all  the  velocity  vectors  lying 
in  a  plane  spanned  by  two  eigenvectors  give  equally  low  errors.  Finally,  if  all  three 
eigenvalues  are  equal,  no  direction  for  f  is  preferred,  since  the  ellipsoid  becomes  the 
unit  sphere. 

The  case  where  exactly  one  eigenvalue  is  zero  also  has  a  simple  geometrical  inter¬ 
pretation.  The  surface  defined  by  g{Uf  Vt  W)  =  0  is  a  straight  line,  which  can  be  seen 
easily  from  the  following  equation: 

X,£/2-fX2K2  =*0 


(47) 
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written  for  the  case  when  X3  is  zero.  (Remember  that  Xi  and  X2  are  both  positive.) 
Clearly  the  unit  sphere  intersects  this  line  in  exactly  two  points,  one  of  them  cor¬ 
responding  to  positive  values  for  Z  in  equation  (34). 

The  method  which  we  just  described  can  be  easily  implemented.  To  this  end,  the 
problem  can  be  discretized.  An  expression  similar  to  (31)  can  be  derived  where  the 
integrals  are  approximated  by  sums.  Our  minimization  method  can  then  be  applied 
to  these  sums.  The  resulting  equations  are  similar  to  ones  described  in  this  section, 
with  summation  replacing  integration.  We  implemented  the  resulting  algorithm  and 
tested  it  using  synthetic  data  including  additive  noise.  Thes  results  agreed  with  our 
expectations. 

One  can  use  the  ratio  of  the  biggest  to  the  smallest  eigenvalue,  the  so  called 
condition  number  [13],  as  a  measure  of  confidence  in  the  computed  velocity.  The  result 
is  very  sensitive  to  errors  in  the  measurements  unless  this  ratio  is  much  bigger  than 
one. 

Curiously,  the  same  error  integral  as  (35)  is  obtained  in  the  case  where  the  MLZuv 
norm  is  used: 


c 


c 


II  f{x,y)  !|z«v=  f  [  [f{x,  y)Z(x, y)]J(ua  -f  v2)dxdy.  (48) 

J—hJ  —w 

We  can  arrive  at  a  similar  solution  by  multiplying  the  integrand  in  (17)  by  Z 2  instead 
of  by  a 7  -f  In  that  case  the  minimization  is  carried  out  with  respect  to  the  MLz 

norm  defined  by: 


/h  /*UI 

/  [f{x,y)Z{x,y)]2dxdy.  (49) 

— til 

Here  optical  flow  velocities  for  points  which  are  further  away  are  weighted  more  heavily. 
This  is  most  appropriate  when  the  measurement  of  larger  velocities  is  less  accurate. 
We  end  up  with  a  quadratic  form  similar  to  g ,  only  the  integrals  for  the  six  constants 
corresponding  to  a,  6,  c,  d,  e,  and  /  are  a  bit  more  complicated.  Curiously  they  only 
depend  on  the  direction  of  the  optical  flow  at  each  point,  not  its  magnitude. 

Also,  other  constraints  could  be  used.  If  we  insist  on  U2-\-V2  =  1,  for  example,  we 
obtain  a  quadratic  instead  of  a  cubic  equation,  and  if  we  use  W  =  1,  a  linear  equation 
only  need  to  be  solved.  The  disadvantage  of  these  approaches  is  that  the  result  is 
sensitive  to  the  orientation  of  the  coordinate  axes.  Clearly,  in  the  case  of  exact  data, 
we  get  the  right  solution  using  any  of  the  constraints  mentioned  above. 

3.4.  Using  a  Different  Constraint 

The  minimization  scheme  discussed  in  the  previous  section  gives  us  a  unique 
solution  in  most  cases  for  the  velocity  vector  T.  Here  we  propose  a  slightly  different 
approach  which  always  gives  us  a  unique  solution.  Note  that  applying  the  first  step  in 
our  minimization  method  gives  us  a  constraint  between  the  values  of  Z,  the  velocity 
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vector  and  the  optical  flow  at  every  point.  We  can  in  addition  aasume  that  Z  =  Zq  ia 
known  at  a  particular  point,  say  (z<>,  Vo)-  Using  the  MLzuv  norm  in  our  scheme,  we 
want  to  minimise: 


f  f  [uZ-(-U  +  zW))8  +  [vZ-(-V  +  yW0lV  +  wJ)^V-  (50) 

J  —  kJ  —  u/ 


Differentiating  (50)  with  respect  to  Z,  and  setting  the  resulting  expression  equal  to  zero, 
we  obtain: 

*  =  =£±2?.  <»t> 

Ua  +  Va 

Thus  we  propose  to  solve  the  minimization  scheme  under  the  following  constraint: 


Zo(ua0  +  t>o)  ~  (“o «  +  *>00)  =  0  (52) 

where  u0  and  Vo  denote  the  components  of  the  optical  flow  measured  at  (zoiVo)-  The 
error  integral  (50)  becomes  after  substituting  (51): 


/:/>- 


VO 


)adxdy 


(53) 


which  is  the  same  as  (35)  and  is  denoted  by  g(U ,  V,  W)  (37).  Thus  we  want  to  minimize: 

g(U,  V,  W)  +  2 \[Zo(u*  +  vl)  -  (u0a  +  vo0)]  (54) 

where  X  denotes  a  Lagrangian  multiplier.  To  determine  U>V  and  W  the  following 
linear  equations  obtained  by  differentiating  (54)  with  respect  to  Ut  V,  W  and  X  have  to 
be  solved: 


> 


■'  ) 


aU  +  dV  +  /W  +  Xu0  =  0 

dU  +  bV  +  tW  +  Xv0  =  0  (55) 

fU  eV  cW  —  X(zoUo  -+•  Vo«o)  “  0 

uqU  +  v0V  —  (zqUo  *+■  =  — Z0(u*  +  wj). 


These  equations  can  be  written  in  the  form: 

G\T\  -  P  (56) 

where  t\  =  (U,  V,  W,  X)T  and  P  =  (0, 0, 0,  —  Zo(u%  +  t>o))T-  Let  ^e  determinant  of 
G\  be  Aot 

A0  =(d8  —  ab)(x0  uo  +  Vo  wo)8  +  (e8  —  bc)ul  +  (/a  —  ac)v\  -f  (57) 

2  ((*  —  bf)uo(x0uQ  +  voWq)  +  (df  —  ae)v0(x0u0  +  Vo^o)  +  («*  —  e/)«oWo]- 


O 
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Assuming  that  Ao  ^  0  we  can  easily  determine  T\  from  (55): 

Tx-Gi'P'.  (58) 

Introducing  the  following  abbreviation: 

^  _  Z0(ug  +  vg) 

“  A  > 

A0 

we  can  give  these  formulae  for  T\: 

U  =  K[uo(bc  —  e2)  +  v0(ef  ccf)  +  +  yoVo)(bf  —  de)] 

V  =  A>0(e/  —  cd)  +  vQ{ac  —  f2)  +  (x0u0  -f  y0vo)(ae  —  df))  (60) 

W  =  K[u0(de  —  bf)  -f  v0{df  —  ae)  +  (^o«o  +  VoVo){d2  —  *&)] 

X  =  IC[ae2  +  d2  +  fc/2  -  ahc  -  2de/]. 

The  disadvantage  of  this  approach  is  that  the  result  depends  upon  the  values  of  the 
optical  flow  at  a  single  point.  To  circumvent  this  problem  we  propose  to  determine 

average  values  for  Ut  V  and  W  in  the  following  manner.  First  note  that  we  are  only 

interested  in  the  ratios  of  U/W  and  V /W  which  obviously  do  not  depend  upon  the 
(unknown)  value  for  Zq.  Equivalently  we  could  determine  the  value  for  K  from  the 
condition  that  T  should  have  unit  length.  Hence  we  can  determine  values  for  U,  V 
and  W  which  depend  only  upon  the  values  of  the  optical  flow  at  a  single  point  and  the 
coefficient  in  the  matrix  G .  If  we  want  to  remove  the  dependence  of  the  result  on  the 
data  at  a  single  point,  we  can  simply  average  the  values  obtained  for  all  image  points. 

In  the  case  where  the  data  is  exact,  i.e.,  where  the  determinant  of  G,  denoted 
b>  detG,  vanishes,  the  solution  for  f  is  the  same  one  as  obtained  using  no  constraint 
in  our  minimizations  scheme.  To  see  this  just  observe  that  in  that  case  X  =  0  as 
X  ==  —K detG).  In  the  case  where  Ao  =  0,  equations  (55)  have  a  solution  only  when 
detG  sbk  0.  We  do  not  have  to  be  concerned  with  the  case  where  Ao  =  0  but  detG  ^  0 
as  we  can  argue  that  equations  (55)  always  have  to  have  a  solution.  Note  that  our 
method  is  based  on  the  condition  that  Z  is  a  certain  function  (51)  of  U,V,  W.  Hence 
(52)  cannot  impose  a  constraint  which  would  be  impossible  to  satisfy. 

4.  Rotational  Case 

Suppose  now  that  the  motion  of  the  camera  is  purely  rotational.  In  order  to 
determine  the  motion  from  optical  flow  we  again  use  a  least-squares  algorithm  with  the 
MLt  norm  described  in  the  previous  section.  Recall  that  in  this  case  the  optical  flow 
is  (see  (8)): 


(59) 


ur  =s  Ax y  —  B{xi  +  1)  -f  Cy  (61) 

vr  =  A(ya  +  1)  —  Bxy  —  Cx. 


d 


16 


We  will  show  now  in  an  analogous  fashion  to  section  3.1  that  two  different  rotations, 
say  u>i  —  and  =  (At,  flj,  C3)T ,  cannot  generate  the  same  optical  flow. 

Assuming  the  converse,  the  following  equations  have  to  hold  for  all  values  of  z  and  y: 

A\xy  —  Bi(xa  -f  1)  +  Ciy  =  A3xy  —  Ba(xa  +  1)  -f-  Cay  (62) 

A\[ya  +  1)  —  B\xy  —  C\x  =  A3{  ya  + 1)  —  Bjxy  —  Cjx 


from  which  we  can  immediately  deduce  that  =  Oa* 

In  general,  the  direction  of  the  optical  flow  at  two  points  and  its  magnitude  at  one 
point  determine  the  motion  of  a  camera  in  pure  rotation  uniquely.  We  choose  instead 
to  minimize  the  following  expression: 


u  —  u,)2  +  (t;  —  vr)2]dxdy. 


(63) 


As  the  motion  is  purely  rotational,  the  optical  flow  does  not  depend  upon  the  distance 
to  the  surface  and  therefore  we  may  omit  the  first  step  in  our  method.  Thus  we  im¬ 
mediately  differentiate  (63)  with  respect  to  A,B  and  C  and  set  the  resulting  expressions 
equal  to  zero: 


[  f  ((u  —  ur)xy  +  (v  —  Vr)(y3  -I-  1  )]dxdy  —  0 
f  I  ((“  ~  «r)(*2  +  1)  +  (t»  —  vT)xy]dxdy  =  0 

J  — III 

/  /  [(u  —  uT)y  —  (v  —  t >r)x]dxdy  =  0. 

J  —hJ — w 


Rewriting  the  above  equations  we  obtain; 

/*»  rw 


/n  pw 

/  [ttxy  -f  v(ya  +  1  )\dxdy  —  I  /  [urzy  +  vT(ya  -f  i)\dxdy 

-hJ — W  J  —  hJ — W 

/h  ami  Ah  rv> 

I  [u(zJ  -f  1)  -f-  vxy]dxdy  =  /  /  (uf (xJ  -|- 1)  -j-  vtxy]dxdy 

-hJ  — HI  J  —  hJ  —w 

/h  AW  ph  pw 

I  (uy  —  vxjdxdy  =  /  /  [ury  —  vrx]dxdy 

-hJ—w  J  —  KJ  — w 

and  expanding  these  equations  yields: 


(64) 


(65) 


) 


) 


SA  +  dB  +  JC 
ZA  +  IB  + 

JA  -j-  2B  ■)■  cC 


(66) 
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where: 


t 


» 


and: 


a  =  f  f  [*  V  +  (V3  +  1  )*]dxdy 

J — hJ—Mf 

5  =  J  +  1)*  +  x*ya]dxdy 

e  =  /  /  (*a  +  V3]d*dy 

J  —kJ—W 

--/.X  [xy(x*  +  ya  -+■  2)]<fxdy 

/■*  r 

Z  =  —  I  I  ydxdy 

?  =  -jXL"“*1 

*  —  [  f  [u*y  +  «(y3  +  l)]dxdy 

J — KJ  — ui 

7  =  —  /  f  [u(x3  +  1)  +  vxy]d*dy 
J - hJ - 10 

fh  =  /  /  [uy  —  vx]dxdy. 

j — — m 


(67) 


(68) 


If  we  call  the  coefficient  matrix  in  (66)  M  and  the  column  vector  on  the  rightrhand  side 
fl,  then  we  have: 

MQ  =  fl.  (69) 

Thus,  provided  the  matrix  Af  is  non~singularf  we  can  compute  the  rotation  as  follows: 

a  =  M-l».  (70) 

It  is  easy  to  see  that  the  matrix  M  is  non-singular  in  the  special  case  of  a  rectangular 
image  plane  since  then: 

a  -  4»Mt  +  —  + 1)  +  — 

b  =»  4toh(—  -f  —  +  1)  H - - — 


Awh 


{w'  +  h3) 

?  =  o. 


(71) 
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So  in  the  case  of  a  rectangular  image  plane,  the  matrix  is  diagonal,  which  makes  it 
particularly  easy  to  compute  its  inverse.  In  fact,  the  matrix  is  diagonal  if  the  image 
plane  is  symmetrical  about  the  x-axis  and  the  y- axis.  As  the  extent  of  the  image  plane 
decreases,  however,  the  matrix  M  becomes  ill  conditioned.  That  is  inaccuracies  in  the 
three  integrals  (£,  7,  and  fh)  computed  from  the  observed  flow  are  greatly  magnified. 
This  makes  sense  since  we  cannot  expect  to  accurately  determine  the  component  of 
rotation  about  the  optical  axis  when  observations  are  confined  to  a  small  cone  of 
directions  about  the  optical  axis.  Again,  in  our  numerical  implementation  of  the 
algorithm  the  integrals  in  (67)  can  be  approximated  by  sums. 


5.  General  Motion 

We  would  like  now  to  apply  a  least- squares  algorithm  to  determine  the  motion  of  a 
camera  from  optical  flow  given  no  a  priori  assumptions  about  the  motion.  It  is  plain  that 
a  least-squares  method  is  easiest  to  use  when  the  resulting  equations  are  linear  in  all  the 
motion  parameters.  Unfortunately,  there  exists  no  norm  which  will  allow  us  to  achieve 
this  goal.  We  did  find  a  norm,  however,  which  resulted  in  equations  that  are  linear  in 
some  of  the  unknowns  and  quadratic  in  the  others.  We  again  attack  the  minimization 
problem  using  the  MLafi  norm  under  the  constraint  that  U2  +  V2  +  W2  =  1.  The 
ensueing  equations  are  polynomials  in  the  unknowns  U ,  V,  Wf  A,  B  and  C  and  can  be 
solved  by  a  standard  iteration  method  like  Newton’s  method  or  Bairstow’s  method  [12] 
or  by  an  interpolation  scheme  like  regul a  falsi  [12].  The  expression  we  wish  to  minimize 
is: 

/J>  -  (f  -  «r)P  +  (V  -  (|  +  «r)]SK«a  +  0*)d*dV'  (72) 

The  first  step  is  to  differentiate  the  integrand  of  (72)  with  respect  to  Z  and  set  the 
resulting  expression  equal  to  aero: 


<*a  +  0a 


(73) 


(u  —  ur)a  +  (v  —  vr)/j‘ 

We  introduce  the  Langrangian  multiplier  X  as  before  and  attempt  to  minimise: 

f  f  [(ti  —  ur)0  —  (v  —  vr)a]adzdy  +  HU3  +  V2  -f  Wa  —  1).  (74) 

J  —v — w 

The  equations  we  have  to  solve  to  determine  the  motion  parameters  are  obtained  by 


o 


;Hv*> 
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differentiation 

rk 


f  f  ((“  —  “r)/J  —  (v  —  vr)a][— xy/3  +  (y2  +  l)a]dxdy  =  0 
f  f  ((«  —  Ur)P  —  {v  —  wr)a][(i3  -I-  1)/?  —  xya]dxdy  —  0 

J  —  hJ — w 

I  I  [(«  —  Ur)P  —  (v  —  vr)a][y/J  +  xajdxdy  =  0 

J  —  hJ — W 

I  I  [(**  —  ur)P  —  («  —  vr)a]vrdxdy  -)-  \U  =  0 

J  —  _t* 

/  /  l(u  —  U*)P  —  (v  —  «r)a]urdidy  +  XV  =  0 

J—kJ—w 

f  f  [(«  —  ur)/J  —  (v  —  vf)a](-«r*  -f  vTy)dxdy  -(-  XIV  =  0 

>  —  KJ  -  w 


u3  +  v*  +  wa  =  1. 


(75) 


Note  that  the  first  three  of  these  equations  are  linear  in  A,  B  and  C  from  which  these 
parameters  can  be  determined  uniquely  in  terms  of  U,  V  and  IV.  Then  we  can  determine 
U,V  and  IV  from  the  last  four  equations  by  a  numerical  method.  To  this  end,  the 
problem  can  be  discretized  and  equations  analogous  to  (75)  derived,  where  summation 
of  the  appropriate  expressions  is  used  instead  of  integration. 


6.  Summary 

Our  objective  was  to  devise  a  method  for  determining  the  motion  of  a  camera  from 
optical  flow  which  allows  for  noise  in  the  measured  data.  The  least- squares  method 
which  we  proposed  in  this  paper  meets  our  goal  and  is  also  suitable  for  numerical 
implementation.  An  important  application  of  our  results  is  in  passive  navigation.  Here 
the  path  and  instantaneous  altitude  of  a  vehicle  is  to  be  determined  from  information 
gleaned  about  the  environment  without  the  emission  of  sampling  radiation  from  the 
vehicle. 
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